﻿using System;

namespace Gmantis.Controls.Chart.Extended
{
    internal class LeastSquaresMethod
    {
        // Methods
        public static double Aproximate(FitType fit, double[] coeffs, double xvalue)
        {
            double num = 0.0;
            switch (fit)
            {
                case FitType.Polynom:
                    {
                        double num2 = 1.0;
                        for (int i = 0; i < coeffs.Length; i++)
                        {
                            if (i > 0)
                            {
                                num2 *= xvalue;
                            }
                            num += coeffs[i] * num2;
                        }
                        return num;
                    }
                case FitType.Exponent:
                    return (coeffs[0] * Math.Exp(xvalue * coeffs[1]));

                case FitType.Logarithmic:
                    return (coeffs[0] * Math.Log(xvalue * coeffs[1]));

                case FitType.Power:
                    return (coeffs[0] * Math.Pow(xvalue, coeffs[1]));

                case FitType.Fourier:
                    num = coeffs[0];
                    for (int j = 1; j < coeffs.Length; j++)
                    {
                        int num5 = (j + 1) / 2;
                        if ((j % 2) == 1)
                        {
                            num += coeffs[j] * Math.Cos(num5 * xvalue);
                        }
                        else
                        {
                            num += coeffs[j] * Math.Sin(num5 * xvalue);
                        }
                    }
                    return num;
            }
            return num;
        }

        private static void ComputeNormalEquations(int np, int nt, double[,] basis, double[] y, double[,] a)
        {
            for (int i = 0; i < nt; i++)
            {
                double num = 0.0;
                for (int j = 0; j < np; j++)
                {
                    num += y[j] * basis[j, i];
                }
                a[i, nt] = num;
                for (int k = i; k < nt; k++)
                {
                    num = 0.0;
                    for (int m = 0; m < np; m++)
                    {
                        num += basis[m, k] * basis[m, i];
                    }
                    a[k, i] = num;
                    a[i, k] = num;
                }
            }
        }

        private static void ComputeStatistics(int np, int nt, double factor, double[] y, double[] solution, double[,] basis, double[] yfit, double[] Residuals, ref double StandardDeviation, ref double Variance, FitType Fit)
        {
            double num = 0.0;
            for (int i = 0; i < np; i++)
            {
                yfit[i] = 0.0;
                for (int j = 0; j < nt; j++)
                {
                    yfit[i] += solution[j] * basis[i, j];
                }
                switch (Fit)
                {
                    case FitType.Exponent:
                        yfit[i] = factor * Math.Exp(yfit[i]);
                        break;

                    case FitType.Power:
                        yfit[i] = factor * Math.Exp(yfit[i]);
                        break;
                }
                Residuals[i] = y[i] - yfit[i];
                num += Residuals[i] * Residuals[i];
            }
            Variance = num;
            if (np == nt)
            {
                StandardDeviation = 0.0;
            }
            else
            {
                StandardDeviation = Math.Sqrt(num / ((double)(np - nt)));
            }
        }

        private static void CreateBasis(int np, double[] x, double[] y, int nt, double[] solution, double[,] basis, FitType fit, bool zero)
        {
            if (np < 2)
            {
                throw new Exception("Incompatible data: Less than 2 data points.");
            }
            if (nt < 1)
            {
                throw new Exception("Incompatible data: Less than 1 coefficient in the fit");
            }
            if (nt > np)
            {
                throw new Exception("Incompatible data: Number of data points less than number of terms");
            }
            for (int i = 0; i < nt; i++)
            {
                solution[i] = 0.0;
            }
            switch (fit)
            {
                case FitType.Polynom:
                    PolyCreateBasisFunctions(np, nt, x, basis, zero);
                    return;

                case FitType.Exponent:
                    ExpoCreateBasisFunctions(np, x, basis);
                    return;

                case FitType.Logarithmic:
                    LogCreateBasisFunctions(np, x, basis);
                    return;

                case FitType.Power:
                    PowerCreateBasisFunctions(np, x, basis);
                    return;

                case FitType.Fourier:
                    FourierCreateBasisFunctions(np, nt, x, basis);
                    return;
            }
        }

        private static void ExpoCreateBasisFunctions(int np, double[] x, double[,] Basis)
        {
            for (int i = 0; i < np; i++)
            {
                Basis[i, 0] = 1.0;
                Basis[i, 1] = x[i];
            }
        }

        private static double ExpoTransform(int np, double[] x, double[] y, double[] xt, double[] yt)
        {
            for (int i = 0; i < np; i++)
            {
                xt[i] = x[i];
            }
            double num2 = (y[0] < 0.0) ? ((double)(-1)) : ((double)1);
            for (int j = 0; j < np; j++)
            {
                double d = num2 * y[j];
                if (d <= 0.0)
                {
                    throw new Exception("Incompatible data: for expo fit all Y values must have the same sign");
                }
                yt[j] = Math.Log(d);
            }
            return num2;
        }

        private static void FourierCreateBasisFunctions(int np, int nt, double[] x, double[,] basis)
        {
            for (int i = 0; i < np; i++)
            {
                basis[i, 0] = 1.0;
                if (nt > 1)
                {
                    basis[i, 1] = Math.Cos(x[i]);
                }
                if (nt > 2)
                {
                    basis[i, 2] = Math.Sin(x[i]);
                }
                for (int j = 3; j < nt; j++)
                {
                    int num3 = (j + 1) / 2;
                    if ((j % 2) == 1)
                    {
                        basis[i, j] = Math.Cos(num3 * x[i]);
                    }
                    else
                    {
                        basis[i, j] = Math.Sin(num3 * x[i]);
                    }
                }
            }
        }

        private static void LogCreateBasisFunctions(int np, double[] x, double[,] basis)
        {
            for (int i = 0; i < np; i++)
            {
                basis[i, 0] = x[i];
                basis[i, 1] = 1.0;
            }
        }

        private static double LogTransform(int np, double[] x, double[] y, double[] xt, double[] yt)
        {
            for (int i = 0; i < np; i++)
            {
                yt[i] = y[i];
            }
            double num2 = (x[0] < 0.0) ? ((double)(-1)) : ((double)1);
            for (int j = 0; j < np; j++)
            {
                double d = num2 * x[j];
                if (d <= 0.0)
                {
                    throw new Exception("Incompatible data: for log fit all Y values must have the same sign");
                }
                xt[j] = Math.Log(d);
            }
            return num2;
        }

        private static void NormalizeAndSolveGauss(int np, int nt, double[,] basis, double[] y, double[] solution)
        {
            double[,] a = new double[nt, nt + 1];
            ComputeNormalEquations(np, nt, basis, y, a);
            if (!SolveGauss(a, solution, 0.0))
            {
                throw new Exception("Incompatible data: No solution");
            }
        }

        private static void PolyCreateBasisFunctions(int np, int nt, double[] x, double[,] basis, bool zero)
        {
            for (int i = 0; i < np; i++)
            {
                if (zero)
                {
                    basis[i, 0] = x[i];
                    if (nt > 1)
                    {
                        basis[i, 1] = x[i] * x[i];
                    }
                }
                else
                {
                    basis[i, 0] = 1.0;
                    basis[i, 1] = x[i];
                }
                for (int j = 2; j < nt; j++)
                {
                    basis[i, j] = x[i] * basis[i, j - 1];
                }
            }
        }

        private static void PowerCreateBasisFunctions(int np, double[] x, double[,] basis)
        {
            for (int i = 0; i < np; i++)
            {
                basis[i, 0] = 1.0;
                basis[i, 1] = x[i];
            }
        }

        private static double PowerTransform(int np, double[] x, double[] y, double[] xt, double[] yt)
        {
            double num = (y[0] < 0.0) ? ((double)(-1)) : ((double)1);
            for (int i = 0; i < np; i++)
            {
                if (x[i] <= 0.0)
                {
                    throw new Exception("Incompatible data: x values should be positive.");
                }
                double d = num * y[i];
                if (d <= 0.0)
                {
                    throw new Exception("Incompatible data: y values must all have the same sign.");
                }
                xt[i] = Math.Log(x[i]);
                yt[i] = Math.Log(d);
            }
            return num;
        }

        private static bool SolveGauss(double[,] a, double[] x, double epsilon)
        {
            int num2;
            int num3;
            int length = a.GetLength(0);
            if ((x.Length < length) || (a.GetLength(1) < (length + 1)))
            {
                throw new ArgumentException("Dimension of matrix is not correct.");
            }
            bool flag = true;
            for (num2 = 0; num2 < length; num2++)
            {
                int num4 = num2;
                double num5 = Math.Abs(a[num2, num2]);
                num3 = num2 + 1;
                while (num3 < length)
                {
                    if (num5 < Math.Abs(a[num3, num2]))
                    {
                        num5 = Math.Abs(a[num3, num2]);
                        num4 = num3;
                    }
                    num3++;
                }
                if (Math.Abs(num5) > epsilon)
                {
                    double num6;
                    num3 = num2;
                    while (num3 <= length)
                    {
                        num6 = a[num2, num3];
                        a[num2, num3] = a[num4, num3];
                        a[num4, num3] = num6;
                        num3++;
                    }
                    for (num4 = num2 + 1; num4 < length; num4++)
                    {
                        num6 = a[num4, num2] / a[num2, num2];
                        a[num4, num2] = 0.0;
                        num3 = num2 + 1;
                        while (num3 <= length)
                        {
                            double num1 = a[num4, num3];
                            num1 -= num6 * a[num2, num3];
                            num3++;
                        }
                    }
                }
                else
                {
                    flag = false;
                    break;
                }
            }
            if (flag)
            {
                for (num2 = length - 1; num2 >= 0; num2--)
                {
                    x[num2] = a[num2, length];
                    for (num3 = num2 + 1; num3 < length; num3++)
                    {
                        x[num2] -= a[num2, num3] * x[num3];
                    }
                    x[num2] /= a[num2, num2];
                }
            }
            return flag;
        }

        public static void SolveLeastSquares(int np, double[] x, double[] y, int nt, double[] solution, FitType fit, bool zero)
        {
            double standardDeviation = 0.0;
            double variance = 0.0;
            SolveLeastSquares(np, x, y, nt, solution, fit, null, null, ref standardDeviation, ref variance, zero, false);
        }

        public static void SolveLeastSquares(int np, double[] x, double[] y, int nt, double[] solution, FitType fit, double[] yfit, double[] residuals, ref double StandardDeviation, ref double Variance, bool zero, bool orto)
        {
            double factor = 0.0;
            double[] xt = x;
            double[] yt = y;
            if (((fit == FitType.Exponent) || (fit == FitType.Logarithmic)) || (fit == FitType.Power))
            {
                nt = 2;
                xt = new double[np];
                yt = new double[np];
                switch (fit)
                {
                    case FitType.Exponent:
                        factor = ExpoTransform(np, x, y, xt, yt);
                        break;

                    case FitType.Logarithmic:
                        factor = LogTransform(np, x, y, xt, yt);
                        break;

                    case FitType.Power:
                        factor = PowerTransform(np, x, y, xt, yt);
                        break;
                }
            }
            if (zero && (fit == FitType.Polynom))
            {
                factor = solution[0];
                nt--;
            }
            double[,] basis = new double[np, nt];
            CreateBasis(np, xt, yt, nt, solution, basis, fit, zero);
            if (orto)
            {
                if (!SolveOrtho(np, nt, basis, yt, solution))
                {
                    throw new Exception("Incompatible data: No solution.");
                }
            }
            else
            {
                NormalizeAndSolveGauss(np, nt, basis, yt, solution);
            }
            if ((yfit != null) && (residuals != null))
            {
                ComputeStatistics(np, nt, factor, y, solution, basis, yfit, residuals, ref StandardDeviation, ref Variance, fit);
            }
            switch (fit)
            {
                case FitType.Polynom:
                    if (zero)
                    {
                        for (int i = solution.Length - 1; i > 0; i--)
                        {
                            solution[i] = solution[i - 1];
                        }
                        solution[0] = factor;
                    }
                    break;

                case FitType.Exponent:
                    solution[0] = factor * Math.Exp(solution[0]);
                    return;

                case FitType.Logarithmic:
                    solution[1] = factor * Math.Exp(solution[1] / solution[0]);
                    return;

                case FitType.Power:
                    solution[0] = factor * Math.Exp(solution[0]);
                    return;

                default:
                    return;
            }
        }

        private static bool SolveOrtho(int n, int m, double[,] a, double[] b, double[] x)
        {
            if ((n <= 0) || (m <= 0))
            {
                return false;
            }
            double[] numArray = new double[n];
            double[,] numArray2 = new double[n, m];
            double[] numArray3 = new double[(m * (m + 1)) / 2];
            for (int i = 0; i < n; i++)
            {
                for (int k = 0; k < m; k++)
                {
                    numArray2[i, k] = a[i, k];
                }
            }
            int index = -1;
            bool flag = true;
            for (int j = 0; j < m; j++)
            {
                double num;
                double num6 = 0.0;
                for (int num7 = 0; num7 < n; num7++)
                {
                    num = numArray[num7] = numArray2[num7, j];
                    num6 += num * num;
                }
                if (num6 == 0.0)
                {
                    flag = false;
                    break;
                }
                numArray3[++index] = num6;
                num = 0.0;
                for (int num8 = 0; num8 < n; num8++)
                {
                    num += numArray[num8] * b[num8];
                }
                x[j] = num;
                for (int num9 = j + 1; num9 < m; num9++)
                {
                    num = 0.0;
                    for (int num10 = 0; num10 < n; num10++)
                    {
                        num += numArray[num10] * numArray2[num10, num9];
                    }
                    numArray3[++index] = num;
                    num /= num6;
                    for (int num11 = 0; num11 < n; num11++)
                    {
                        double num1 = numArray2[num11, num9];
                        num1 -= numArray[num11] * num;
                    }
                }
            }
            if (flag)
            {
                int num12 = m + 1;
                int num13 = m;
                while (--num13 >= 0)
                {
                    int num14 = index - num13;
                    double num15 = x[num13];
                    for (int num16 = num13 + 1; num16 < m; num16++)
                    {
                        num15 -= numArray3[num16 + num14] * x[num16];
                    }
                    x[num13] = num15 / numArray3[index];
                    index += num13 - num12;
                }
            }
            return flag;
        }
    }
}
